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Jamming transition in a two-dimensional open granular pile with rolling 

resistance 

C. F. M. Magalhaes,^ A. P. F. Atman,G. Combe,J. G. MoreiraP 


We present a molecular dynamics study of the jamming/unjamming transition in two- 
dimensional granular piles with open boundaries. The grains are modeled by viscoelastic 
forces, Coulomb friction and resistance to rolling. Two models for the rolling resistance 
interaction were assessed: one considers a constant rolling friction coefficient, and the 
other one a strain dependent coefficient. The piles are grown on a finite size substrate and 
subsequently discharged through an orifice opened at the center of the substrate. Varying 
the orifice width and taking the final height of the pile after the discharge as the order 
parameter, one can devise a transition from a jammed regime (when the grain flux is always 
clogged by an arch) to a catastrophic regime, in which the pile is completely destroyed 
by an avalanche as large as the system size. A finite size analysis shows that there is a 
finite orifice width associated with the threshold for the unjamming transition, no matter 
the model used for the microscopic interactions. As expected, the value of this threshold 
width increases when rolling resistance is considered, and it depends on the model used 
for the rolling friction. 


I. Introduction 

Granular materials are ubiquitous either in nature 
—desert dunes, beach sand, soil, etc.— or in indus- 
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trial processes as mineral extraction and process¬ 
ing, or food, construction and pharmaceutical in¬ 
dustries m- In fact, any particulate matter made 
of macroscopic solid elements can be classified as 
granular material. The vast phenomenology ex¬ 
hibited by these systems combined with an incom¬ 
plete understanding about the microscopic physical 
mechanisms responsible for the macroscopic behav¬ 
ior of these materials have motivated the increasing 
interest of the physics community in the past years 
Il[5]. 

Although materials of this class are not sensitive 
to thermal fluctuations, they can be found at gas, 
liquid or solid phases [B]. The transition between 
solid and liquid phases in granular matter, which 
is commonly referred to as jamming/unjamming 
transition, has been extensively studied from both 
theoretical and experimental perspectives gUMII]. 
Currently, a great effort is being made to under- 
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stand the nature of this transition, which is still a 
subject of debate m- The jamming/unjamming 
transition is not a specific property of granular 
matter, being observed in many kinds of materi¬ 
als, such as foams [13], emulsions [14], colloids [15], 
gels [12], and also in usual molecular liquids 
—glass transition. Liu and Nagel [5] proposed a 
general phase diagram as an attempt to unify the 
several approaches to study jamming/unjamming 
in disordered materials. This work has motivated 
several theoretical, experimental and numerical in¬ 
vestigations, but a comprehensive understanding of 
this transition is still lacking. 

O’Hern et al. [nm] have performed numeri¬ 
cal simulations of granular materials approaching 
to jamming in two and three dimensions. They 
have explored the packing fraction axis of the gen¬ 
eral phase diagram proposed by Liu and Nagel and 
have demonstrated, by means of finite-size analy¬ 
sis, the existence of a unique critical point in which 
the system jams in the thermodynamic limit. The 
authors have also shown some evidence that this 
point is an ordinary critical point, indicating that 
the jamming transition would be a second order 
phase transition. These results were corroborated 
later by experiments (c.f. Majmudar et al. |10j 'l 
and by simulations (c.f. Manna and Khakhar [20] 1. 
In Ref. [10], the authors have revealed evidence of 
self-organized criticality (SOC) by measuring the 
internal avalanches resulting from the opening of 
an orifice at the bottom of granular piles. 

Experimental investigation of the jamming tran¬ 
sition in granular materials under gravitational 
field has been conducted in a variety of ways, ad¬ 
dressing the role of many parameters, like the grain 
shape, the friction coefficient, and the system ge¬ 
ometry. However, a common feature of all these ap¬ 
proaches is the analysis of the granular flow through 
bottlenecks. 

The jamming of three-dimensional piles seems to 
be settled after the work of Zuriguel et al. m- 
They have demonstrated experimentally, for piles 
composed of different kinds of grains, the diver¬ 
gence on the mean internal avalanche size. It means 
that, as the outlet size approaches a critical value, 
the internal avalanche increases without limit and 
a permanent flow is established. This critical outlet 
is insensitive to the density, stiffness and roughness 
of the grains, but shows a significant dependence on 
the grain shape. For spherical grains, a critical out¬ 


let width Wc ~ 4.94d was obtained, for cylindrical 
grains Wc S.OSd and for rice grains Wc ^ 6.15d. 

Nevertheless, the jamming transition in two- 
dimensional piles is still a question under debate. 
In order to address it, To et al. [22] have carried out 
experiments using two-dimensional hoppers in or¬ 
der to find a critical outlet size for jamming events. 
The jamming probability J has presented a rapid 
decay from J = I to J = 0 close to the aperture 
width w ~ 3.8d, signaling a possible phase transi¬ 
tion. The authors discussed, based on a restricted 
random walk model, the connection between jam¬ 
ming and the arch formation mechanism. Never¬ 
theless, the point needs further investigation, espe¬ 
cially at the limit of high hopper angle. There exist 
several works [231 - 125] focused on the mechanisms of 
arch formation, but none had explored its relation 
to jamming probability. 

Janda et al. [22] have made some progress 
by simulating discharges of two-dimensional silos. 
They have improved the definition of jamming so 
that the internal avalanche size is considered. This 
modification addresses the extremely long relax¬ 
ation times associated with jamming. Within the 
framework of a probabilistic theory concerning the 
arch formation, the authors have tested two hy¬ 
pothesis for the internal avalanche behavior: a 
functional form that predicts a divergence in the 
mean size, and a functional form where the mean 
size exists for all values of the orifice width w. Since 
the latter one is more compatible with the arch for¬ 
mation model, the authors claimed that “no critical 
opening size exists beyond which there is not jam¬ 
ming” [22] • 

The results mentioned so far are related to fully 
confined systems. Recently, a simulation study 
on discharges of granular piles with open bound¬ 
aries m provided new insights on the problem. 
The piles are composed by homogeneous disks in¬ 
teracting via elastic and frictional forces. Using 
finite size analysis, the authors have shown that 
a catastrophic regime, in which the pile is com¬ 
pletely destroyed by the opening of the orifice, is 
well defined. At the limit of infinitely large sys¬ 
tems, the catastrophic regime coincides with the 
unjammed phase, since it implies a divergent in¬ 
ternal avalanche. Hence, the results indicate the 
existence of the jamming transition. It is impor¬ 
tant to note, however, that the pile geometry could 
probably play a role in the causes for this distinct 
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behavior, due to the absence of the Janssen effect, 
but further investigations are necessary to confirm 
this point. 

In the present work, the investigation of jamming 
in 2D open systems is extended in order to con¬ 
sider a rolling resistance term in the grain interac¬ 
tion model, following the prescription adopted by 
Chevalier et al. [55] • The main objective of this 
study is the verification of rolling resistance influ¬ 
ence on jamming. Many factors contribute to the 
appearance of rolling friction, including microslid- 
ing, plastic deformation, surface adhesion, grain 
shape, etc., but mainly, it is due to the contact de¬ 
formation [55] . Here, it will be taken into account 
only the effect due to the contact deformation by 
implementing the micromechanical model proposed 
by Jiang et al. [30]. The rolling friction produces 
a resistance to roll which, among other effects, is 
responsible for granting more stability to granular 
piles ED, and for the occurrence of different types 
of failure modes in granular matter mm- Rolling 
friction was also used to model a system of polygo¬ 
nal grains by making a correspondence between the 
rolling stiffness and the number of sides of the poly¬ 
gon ED . It was demonstrated that it is an essential 
ingredient to reproduce experimental compression 
tests in mixtures of two-dimensional circular and 
rectangular grains ESI- These facts suggest that 
jamming could be affected by rolling friction. Nev¬ 
ertheless, most studies on granular materials based 
on computer simulations do not deal with it. 

The paper is organized as follows: after a review 
in the Introduction, the next section is concerned 
with the methodology. Then, we present the re¬ 
sults and a brief discussion. Finally, the last section 
gathers the conclusions and some perspectives. 

II. Methods 

The jamming transition is assessed by means of 
discharges of granular piles, simulated using the 
molecular dynamics method |36j with the Velocity- 
Verlet algorithm m- In a few words, the molecu¬ 
lar dynamics consists in integrating numerically the 
equations of motion that governs the system dy¬ 
namics. The system is constituted by N free grains 
governed by Newton’s second law and by a finite 
and horizontally aligned substrate made of fixed 
grains. The free grains, which will form the pile, are 


homogeneous bi-dimensional disks that are free to 
translate and rotate around their center, and whose 
radii are uniformly distributed around an average 
value d, a small polydispersity of 5% was imposed 
in order to avoid crystallization effects. All spatial 
quantities will be expressed in terms of d. Since the 
grains have all the same mass density, their masses 
are proportional to their respective areas. Normal¬ 
ized by the mass of the heaviest grain, the masses 
are given by rrii = where dmax stands for 

the diameter of the largest grain. The finite sub¬ 
strate of length L is composed by fixed grains of 
the same kind but with a smaller and fixed diame¬ 
ter ds = O.Id. These grains are aligned horizontally 
and are equally spaced in order to form a grid with¬ 
out gaps. 

The grains are subject to a uniform gravita¬ 
tional field orthogonal to the substrate line, and to 
short range binary interactions. Besides the visco¬ 
elastic and coulomb friction interactions used in 
past works 155115711551 1^5]. two grains in contact are 
also subject to a rolling resistance moment due to 
the finite contact length 1^. The rolling resistance is 
introduced through a micromechanical model of the 
contact line between two grains ES]- This model 
treats the contact as an object formed by a set of 
springs and dashpots connecting the borders of the 
two grains. As one grain rolls over the other, the 
springs in one side of the contact line contract while 
the springs in the opposite side stretch. This config¬ 
uration generates an unbalanced force distribution 
and a consequent moment with respect to the grain 
center (see Fig. [IJ. This moment grows linearly 
as the grain rolls, until the rolling displacement 5r 
reaches some threshold value at which the springs 
located near one end of the contact line break up 
and new ones emerge at the other end. At that 
time, the moment saturates at some value that de¬ 
pends on the properties of the grain. The rolling 
displacement is defined by 5r = “ ujj)At, 

where uJi refers to the angular velocity of grain i, 
At is molecular dynamics time step, and the sum 
runs over time during the whole existence of con¬ 
tact. Based on these assumptions, the authors have 
derived an analytic expression for the rolling resis¬ 
tance moment as a function of rolling displacement. 
They have also proposed a simplified version - the 
one used in this study - in order to improve numer¬ 
ical computations: 
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kr^r-) ^ f^rfel 

_ Sr „ i- lA I ^ f 

\dr\ 

where kr is the rolling stiffness, Hr is the coefficient 
of rolling resistance, and fei is the compressive elas¬ 
tic force, normal to contact line. As can be noted, 
the expression for the rolling resistance momentum 
possesses a striking resemblance with the Coulomb 
static friction force, and as a matter of fact, the 
rolling resistance interaction is implemented in the 
molecular dynamics algorithm in the same way as 
the static friction. The model predicts that the 
rolling stiffness and the coefficient of rolling resis¬ 
tance are related to the contact length Ic by the 
equations kr = knl^ and ^r = where kn and 
fj, are respectively the normal elastic constant and 
friction coefficient. The values used for these pa¬ 
rameters were the same as in Ref. m, fj, = 0.5 and 
kn = 1000 in normalized unities (see [40] for further 
details). Two cases were investigated: systems in 
which the rolling resistance parameters kr and fj,r 
vary according to the above-mentioned equations, 
and systems with fixed rolling resistance parame¬ 
ters, assuming that all contacts have the same de¬ 
formation value. 

The simulation procedure consists of two steps: 
(1) the formation of the granular pile with open 
boundaries, by deposing grains from rest, under 
gravity, over the substrate until an stationary state 
is reached; (2) the discharge itself, which consists in 
opening an orifice of a given width at the center of 
the substrate. In the first step, the initial positions 
of the free grains are randomly sorted along a hori¬ 
zontal line located at height L from the substrate - 
the releasing height is equal to the substrate length. 
To avoid initial overlapping of grains, a 50% filling 
ratio was imposed to each line of grains released, 
and the time interval between successive rows is 
the inverse of the frequency /. Each row was re¬ 
leased after the predecessor had fallen a distance 
equivalent to the maximum grain diameter. This 
deposition protocol mimics a dense rain of grains. 
During deposition, the grains may leave the sys¬ 
tem through the lateral boundaries, so that the to¬ 
tal number of grains in the pile fluctuates as the 
process evolves. The release of grains ceases when 
the number of grains in the pile reaches a station¬ 
ary value. The deposition phase ends only when a 



Figure 1: In panel (a), a perfectly rigid grain rolling 
over another one is exhibited. The contact force is 
concentrated on one point and is aligned through 
the grain center, thus, not generating momentum 
with respect to the center. Panel (b) shows the 
same situation but with deformable grains. Now, 
the contact force is non-uniformly distributed over 
a segment of length Ic (contact length). The re¬ 
sultant force A.c is dislocated from the center line 
and an opposing moment with respect to the center 
appears. 

mechanical equilibrium state is attained, and the 
configuration is recorded for later analysis. The 
equilibrium state must satisfy the following crite¬ 
ria: mechanical stability, absence of slipping con¬ 
tacts, vertical and horizontal force balance, and 
vanishingly small kinetic energy |43] . In the sec¬ 
ond step, the configuration recorded is loaded and 
an outlet of width w is opened at the center of the 
substrate, allowing the grains to flow through it. 
As the grains pass through the outlet, they are re¬ 
moved from the system, in order to improve compu¬ 
tational efficiency. The simulation runs until a new 
equilibrium state is reached, which happens either 
due to the formation of an arch above the outlet or 
after the pile has been completely discharged. The 
remaining pile configuration is then recorded. 

The average height h of the resulting pile after 
discharge, measured from the center of the sub¬ 
strate, was taken as an order parameter to distin¬ 
guish between the two regimes. If h does not change 
significantly with respect to the original height, it 
means that an arch does readily clog the flux af¬ 
ter the orifice opening, but if /i = 0, it means that 
the pile has collapsed, and a jammed state was not 
attained for the corresponding orifice diameter and 
substrate length. As the orifice width w approaches 
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Figure 2: Diagram of a possible equilibrium con¬ 
figuration if rolling resistance is considered. The 
grain on top has only one contact and, even so, it 
can sustain a stable position. The maximum value 
of the angle 9 depends on the friction coefficient 
and on the rolling resistance parameters. 

a certain threshold value Wt, the system suffers a 
transition from jammed to unjammed state. This 
threshold is defined as the orifice width for which 
the height fluctuations is maximum. The connec¬ 
tion to the jamming transition occurs when the 
substrate length diverges since the collapse of an 
infinite substrate pile implies a continuous flowing 
state. Then, the jamming transition can be charac¬ 
terized by a critical aperture width, as defined by 
the following expression 

Wc = lim Wt{L) . (2) 

L—¥00 

III. Results and Discussion 

Two different models of rolling resistance interac¬ 
tion were tested in the simulations: the original 
model proposed by Jiang et al. [30] mentioned 
earlier, with a rolling constant and a coefficient 
of rolling resistance that depends on the contact 
length {kr = knl^ and Mr = k'lc) and a simpler 
derived version, in which these parameters are con¬ 
stant over time assuming that Ic = 0.05 d for all 
contacts. This fixed Ic model assumes a mean con¬ 
tact length equivalent to 5% of the average grain 
diameter, a scenario which could be associated to a 
system composed by polygonal grains, for example, 
a extreme rolling resistance regime. 

For either models, the numerical simulations de¬ 
scribed in the last section were carried out for 



Figure 3: Images of the pile equilibrium states for 
the three grain interaction models before the dis¬ 
charge step. The piles were grown over a L = lOOd 
substrate and are representative samples from each 
type of pile. As shown in the figure, the top image 
represents a pile composed by grains with a fixed 
kr rolling resistance interaction, the image in the 
middle is a pile of grains with a kr 1^ rolling re¬ 
sistance interaction, and the bottom image is a pile 
of grains without rolling resistance. 

various system sizes and orifice widths. Figure 
[3] shows equilibrium configurations of typical piles 
with L = lOOd for the two tested rolling resistance 
models and for the absence of rolling resistance 
case. It can be seen that the inclusion of the rolling 
resistance term modifies significantly the macro¬ 
scopic features of the pile. It provides more stabil¬ 
ity to the structure, which is reflected by a steeper 
free surface, a fact also observed elsewhere |3T]. In¬ 
deed, it should be noted that the rolling resistance 
makes possible some otherwise very unstable two- 
grain configuration, as exemplified in Fig. [2] 

The behavior of the order parameter h as func¬ 
tion of w is presented in Fig. 2] for the three cases. 
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Figure 4: Order parameter h as a function of the 
orifice width w for all types of pile. The graphs were 
generated from the simulation data of L = lOOd 
piles. While the symbols represent the parameter 
h itself, the lines indicate the corresponding fluc¬ 
tuations. The thick line is related to the fixed kr 
curve, the medium thickness line to the kr ~ lc2 
curve, and the dashed line to kr = 0 curve. 

Note that the transition region translates to the 
right as the rolling resistance becomes more impor¬ 
tant, which is an expected result since the increas¬ 
ing of stability allows the formation of larger arches. 
Figure [5] exhibits the dependence of Wt (fluctuation 
maximum) on the system size for the three mod¬ 
els considered and, again, the curves are dislocated 
along the Wt axis. Nevertheless, they all share the 
same functional aspect, which is an evidence that 
the rolling resistance does not change the nature 
of the jamming transition, only the critical value 
of the threshold width. Despite the tendency to a 
heavy tail distribution observed for large L and w, 
in all scenarios, the data suggest the existence of 
the jamming transition. It means that the features 
described in Ref. m to characterize the transi¬ 
tion were also observed in both scenarios tested 
for the rolling resistance. The fitting values were 
Wc = (5.3 ± 0.1)c? for contact length dependent 
kr model and Wc = (7.6 ± 0.2)d for the fixed kr 
one, while for the model without rolling resistance, 
Wc = (5.0 ± 0.1)d [57]. These results indicate that 
the rolling resistance only slightly changes the crit¬ 
ical width when included in a system of disks, ex¬ 
pressed by the contact length model. But, in other 



Figure 5: Graphs of the threshold orifice width wt 
as a function of the reciprocal system size 1/L for 
models with and without rolling resistance. As the 
caption indicates, the symbols represent the data 
obtained from numerical simulations and the lines 
are fitting curves, which provide the respective val¬ 
ues of Wc- 

approaches, as in the fixed length model, it can al¬ 
ter significantly the critical aperture width. 

The probability density functions (PDF) of h for 
the absent rolling friction and for the fixed rolling 
friction models are exhibited in Fig. |6| For each 
model, the PDFs are obtained for three characteris¬ 
tic values of w: w = l.Od , w ^ Wt, and w > Wt- It 
can be noted that the PDFs have an approximately 
Gaussian peak around the initial height for all val¬ 
ues of w, meaning that there is always a certain 
amount of samples that are simply not disturbed 
by the outlet. Apart from these samples, the great 
majority is completely discharged for w > Wt- This 
result is in consonance with that obtained earlier 
in a different context m- In that study, it was 
shown that for large w, all blocked events occurred 
after only few grains had passed through the outlet. 
These facts indicate that the initial conditions may 
play an important role in jamming experiments. 
However, this issue needs further investigation. 

IV. Conclusions 

Evidence of the jamming transition was observed 
in molecular dynamics simulations of open granu¬ 
lar piles with rolling resistance, in consonance with 
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Figure 6: Probability density functions of the order 
parameter h for different values of w in the case of 
absent rolling friction grains (top) and fixed rolling 
friction grains (bottom). 

the results found in granular piles without this in¬ 
teraction. This result strengthens the expectation 
that the transition exists in real granular piles and 
is probably affected by the system geometry. We 
observe that when there was rolling resistance, the 
piles built were more stable, denoted for the large 
mean height of the samples, and also more ro¬ 
bust against perturbations, since the critical aper¬ 
ture width increased. In future works, we plan to 
present a detailed study of the arching statistics for 
the two approaches considered here. 
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